% Reading gridded weather data and matching it with lat/lon from ATUS.
% To run this code on future years, you will need to edit the year list and
% the list of leapyears.

%% Preliminaries
close all;
clear all;
temp_dir = '/Volumes/Data 2TB/data/weather';
cd '/Users/jgs/Dropbox/research/projects/active/ATUS_project/Data';

%% For each year, read in the weather data, match it with the ATUS lat/lon
years = 2003:2013;
atus_round = csvread('atus_lat_lon_for_weather.csv');

% Get ready to store the daily temperature data
% Note that 2004, 2008, and 2012 were leap years
atus_temp = NaN(length(atus_round(:,1)),365*length(years)+3);
leap_year = [365,366,365,365,365,366,365,365,365,366,365,365] - 365;

% Loop over temperature datasets
start_day = 1;
end_day = 365;
for y = 1:length(years);
    % Read data
    read_netcdf(strcat(temp_dir, '/air.sig995.', int2str(years(y)),'.nc'));
    % Match with given ATUS location

    for i = 1:length(atus_round(:,1));
        atus_temp(i,start_day:end_day) = air(lat==atus_round(i,1),lon==atus_round(i,2),:);
    end;
    fprintf('%4.0f\n',year(y));
    start_day = end_day + 1;
    end_day = end_day + 365 + leap_year(y+1);
end;

% Reattach lat-lon
atus_out = [atus_round, atus_temp];
% Output file
csvwrite('atus_temperature.csv',atus_out);




%% Just to check things, make a map of a single day
% To make this map, load at least one 
R_all = spatialref.GeoRasterReference; 
R_all.RasterSize = [73 144]; 
R_all.Latlim = [-90 90]; 
R_all.Lonlim = [0 357.5]; 
%R_all.ColumnsStartFrom = 'north';

h = figure;
axesm('MapProjection','eckert4','MapLonLimit',[0 360]);
geoshow(air(:,:,139), R_all,'DisplayType', 'texturemap');
hcb = colorbar('southoutside');

% Plotting ATUS data
%h1 = figure;
%axesm('MapProjection','eckert4','MapLonLimit',[0 360]);
%geoshow(atus_data,'DisplayType', 'texturemap');
